
#################################################################################
# File Name:  	appendix.R 
# Paper:			  Beliefs about Climate Action Consequences under Weak Global Institutions: #               Sectors, Home Bias, and International Embeddedness
# Author: 		  Patrick Bayer (Strathclyde) and Federica Genovese (Essex)
#
# Purpose:      Main results
# Data input:   "./gep_data.csv"
# Date:         27 September 2020
#
#################################################################################

rm(list=ls())

library(ggplot2)
library(dplyr)
library(ggpubr)
library(reshape2)


# Load data
f<-read.csv("gep_data.csv")


# Drop incomplete responses
f <- f[f$Progress==100,]

# Exclude those who failed the attention check
f <- f[f$attentioncheck==3 & is.na(f$attentioncheck)==FALSE,]

# Exclude pacers who were faster than 3 minutes 
f<- f[f$Durationinseconds>180,]
mean(f$Durationinseconds)/60

# Clean data (Don't know responses)
f$devep_ccperformance[f$devep_ccperformance==5] <- NA
f$deving_ccperformance[f$deving_ccperformance==5] <- NA
#f$Commonwealth[f$Commonwealth==4] <- NA

# ================================================================
# Variable creation
# ================================================================

# Create indicators for experimental conditions
f$treat_vign[is.na(f$t_control_Page.Submit)==FALSE] <- 0
f$treat_vign[is.na(f$t_t1_domwin_Page.Submi)==FALSE] <- 1
f$treat_vign[is.na(f$t_t2_domlos_Page.Submit)==FALSE] <- 2
f$treat_vign[is.na(f$t_t3_forwin_Page.Submit)==FALSE] <- 3
f$treat_vign[is.na(f$t_t4_forlos_Page.Submit)==FALSE] <- 4

f$treat_vign <- factor(f$treat_vign, levels=c(0,1,2,3,4), labels=c("C", "Tr_DomWin", "Tr_DomLose", "Tr_ForWin", "Tr_ForLose"))

table(f$treat_vign, useNA="ifany")


# ================================================================
# Appendix A1: Descriptive Statistics
# ================================================================

library(reporttools)

summary(f)

# Define descriptive stats variables
desc.demo <- f[, c("gender","age","residence","gor_england", "education", "income", "empl_sector")]
desc.pol <- f[, c("political_id", "X2017_election", "X2017_election_party", "parties", "Brexit", "Commonwealth", "travels")]
desc.climate <- f[, c("cc_seriousproblem", "cc_concern", "UK_ccperformance", "devep_ccperformance", "deving_ccperformance", "heardparisaccord")]
desc.sector <- f[, c("sectorperformance", "transportsectorfeel", "energysectorfeel")]
desc.outcomes <- f[, c("vign_sectoralwin", "vign_beliefUKwins", "vign_beliefGerwins", "vign_beliefdevingwin", "vign_supportUKaction", "vign_supportdevelopd", "vign_supportdeveling")]

# Produce descriptive stats tables
tableNominal(vars=desc.demo, cap="Descriptive statistics", lab="tab:desc.demo", longtable=FALSE)
tableNominal(vars=desc.pol, cap="Descriptive statistics", lab="tab:desc.pol", longtable=FALSE)

tableContinuous(vars=desc.sector, cap="Descriptive statistics", lab="tab:desc.sector", longtable=FALSE, stats=c("n", "min", "median", "mean", "max", "na"))
tableContinuous(vars=desc.outcomes, cap="Descriptive statistics", lab="tab:desc.outcomes", longtable=FALSE, stats=c("n", "min", "median", "mean", "max", "na"))


# ================================================================
# Appendix A2: Balance Table
# ================================================================


# Call variable from descriptive stats one-by-one
var <- colnames(desc.climate)
OUT <- list()

for (j in 1:length(var)){

  pvallist <- NA
  for (i in 2:5) {
    temp.data <- f[as.numeric(f$treat_vign)==1 | as.numeric(f$treat_vign)==i,]
    pvallist[i-1] <- t.test(temp.data[,colnames(temp.data)==var[j]] ~ treat_vign, data=temp.data)$p.value
  }
  OUT[[j]] <- c(var[j],p.adjust(pvallist,method="BH"))
}

full = do.call(rbind, OUT)
full

# ================================================================
# Appendix A3: Main Results Table
# ================================================================

# Please see 'main.R' file 


# ================================================================
# Appendix A4: Support for Climate Action
# ================================================================

# Support for UK climate action
h1 <- group_by(f, treat_vign) %>%
  summarise(
    count = n(),
    mean = mean(vign_supportUKaction, na.rm = TRUE),
    sd = sd(vign_supportUKaction, na.rm = TRUE)
  )


# ANOVA: UK climate action
res.aov <- aov(vign_supportUKaction ~ treat_vign, data = f)
summary(res.aov)
TukeyHSD(res.aov)

# Multiple comparison p-value correction
pvallist <- NA
for (i in 2:5) {
  pvallist[i-1] <- t.test(vign_supportUKaction ~ treat_vign, data=f[as.numeric(f$treat_vign)==1 | as.numeric(f$treat_vign)==i,])$p.value
}
p.adjust(pvallist,method="BH")


# Support for German climate action
h2 <- group_by(f, treat_vign) %>%
  summarise(
    count = n(),
    mean = mean(vign_supportdevelopd, na.rm = TRUE),
    sd = sd(vign_supportdevelopd, na.rm = TRUE)
  )


# ANOVA: Germany climate action
res.aov <- aov(vign_supportdevelopd ~ treat_vign, data = f)
summary(res.aov)
TukeyHSD(res.aov)

# Multiple comparison p-value correction
pvallist <- NA
for (i in 2:5) {
  pvallist[i-1] <- t.test(vign_supportdevelopd ~ treat_vign, data=f[as.numeric(f$treat_vign)==1 | as.numeric(f$treat_vign)==i,])$p.value
}
p.adjust(pvallist,method="BH")


# ================================================================
# Figure A1: Support for UK/German climate action
# ================================================================


# Restructure data for better plotting options
df.plot <- cbind(f[,c("treat_vign","vign_supportUKaction","vign_supportdevelopd")])
colnames(df.plot)[2:3] <- c("Support: UK climate action", "Support: German climate action")
g <- melt(df.plot,value.name="outcome")   

pdf(file="./support_UKGER.pdf", width=8, height=5)

cols <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")

# Construct data set for means
meandf <- data.frame(variable=c("Support: UK climate action", "Support: German climate action"), mean=c(h1$mean[1],h2$mean[1]))

ggerrorplot(g, x="treat_vign", y="outcome", color="variable",
            desc_stat="mean_ci", 
            error.plot="errorbar",
            add="mean",
            legend="none",
            ylab="DV (ordered)", xlab = "Experimental conditions",
            palette=c(cols[2],cols[3]),
            ylim=c(2,4)
            ) +
            facet_grid(~variable, scales="free") +
            scale_x_discrete(breaks=c("C", "Tr_DomWin", "Tr_DomLose", "Tr_ForWin", "Tr_ForLose"), labels=c("Control", "T1: UK \n sector wins", "T2: UK \n sector loses", "T3: German \n sector wins", "T4: German \n sector loses")) +
            geom_hline(data=meandf, aes(yintercept=mean), linetype=c("dashed","dashed"), size=.5,col=c(cols[2],cols[3])) +
            theme_bw() +
            theme(text = element_text(size=11),
            axis.text.x = element_text(size=9),
            legend.position="none"
            )
dev.off()



# ================================================================
# Appendix A5: Analysis for Developing Countries
# ================================================================


# Belief that UK wins
h1 <- group_by(f, treat_vign) %>%
  summarise(
    count = n(),
    mean = mean(vign_beliefUKwins, na.rm = TRUE),
    sd = sd(vign_beliefUKwins, na.rm = TRUE)
  )

# Belief that Germany wins
h2 <- group_by(f, treat_vign) %>%
  summarise(
    count = n(),
    mean = mean(vign_beliefGerwins, na.rm = TRUE),
    sd = sd(vign_beliefGerwins, na.rm = TRUE)
  )


# Belief that India wins
h3 <- group_by(f, treat_vign) %>%
  summarise(
    count = n(),
    mean = mean(vign_beliefdevingwin, na.rm = TRUE),
    sd = sd(vign_beliefdevingwins, na.rm = TRUE)
  )


# ANOVA: India wins
res.aov <- aov(vign_beliefdevingwin~ treat_vign, data = f)
summary(res.aov)
TukeyHSD(res.aov)


# Multiple comparison p-value correction
pvallist <- NA
for (i in 2:5) {
  pvallist[i-1] <- t.test(vign_beliefdevingwin ~ treat_vign, data=f[as.numeric(f$treat_vign)==1 | as.numeric(f$treat_vign)==i,])$p.value
}
p.adjust(pvallist,method="BH")


# ================================================================
# Figure A2: Climate beliefs for the UK, Germany, and India
# ================================================================

# Restructure data for better plotting options
df.plot <- cbind(f[,c("treat_vign","vign_beliefUKwins","vign_beliefGerwins", "vign_beliefdevingwin")])
colnames(df.plot)[2:4] <- c("Belief: UK wins", "Belief: Germany wins", "Belief: India wins")
g <- melt(df.plot,value.name="outcome")                     


pdf(file="./belief_all_win.pdf", width=8, height=8)

cols <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")

# Construct data set for means
meandf <- data.frame(variable=c("Belief: UK wins", "Belief: Germany wins", "Belief: India wins"), mean=c(h1$mean[1],h2$mean[1], h3$mean[1]))


ggerrorplot(g, x="treat_vign", y="outcome", color="variable",
            desc_stat="mean_ci", 
            error.plot="errorbar",
            add="mean",
            legend="none",
            ylab="DV (ordered)", xlab = "Experimental conditions",
            palette=c(cols[2],cols[3],cols[4]),
            ylim=c(1.5,3.5)
            ) +
            facet_wrap(~variable, scales="free", ncol=2) +
            scale_x_discrete(breaks=c("C", "Tr_DomWin", "Tr_DomLose", "Tr_ForWin", "Tr_ForLose"), labels=c("Control", "T1: UK \n sector wins", "T2: UK \n sector loses", "T3: German \n sector wins", "T4: German \n sector loses")) +
            geom_hline(data=meandf, aes(yintercept=mean), linetype=c("dashed","dashed","dashed"), size=.5,col=c(cols[2],cols[3],cols[4])) +
            theme_bw() +
            theme(text = element_text(size=11),
              axis.text.x = element_text(size=9),
            legend.position="none"
            )
dev.off()



# Support for Indian climate action
group_by(f, treat_vign) %>%
  summarise(
    count = n(),
    mean = mean(vign_supportdeveling, na.rm = TRUE),
    sd = sd(vign_supportdeveling, na.rm = TRUE)
  )


# ANOVA: UK climate action
res.aov <- aov(vign_supportdeveling ~ treat_vign, data = f)
summary(res.aov)
TukeyHSD(res.aov)

# Multiple comparison p-value correction
pvallist <- NA
for (i in 2:5) {
  pvallist[i-1] <- t.test(vign_supportdeveling ~ treat_vign, data=f[as.numeric(f$treat_vign)==1 | as.numeric(f$treat_vign)==i,])$p.value
}
p.adjust(pvallist,method="BH")

# ================================================================
# Sectoral beliefs (footnote in main text)
# ================================================================


# Belief that sector wins
group_by(f, treat_vign) %>%
  summarise(
    count = n(),
    mean = mean(vign_sectoralwin, na.rm = TRUE),
    sd = sd(vign_sectoralwin, na.rm = TRUE)
  )

# ANOVA: UK wins
res.aov <- aov(vign_sectoralwin ~ treat_vign, data = f)
summary(res.aov)
TukeyHSD(res.aov)

# Multiple comparison p-value correction
pvallist <- NA
for (i in 2:5) {
  pvallist[i-1] <- t.test(vign_sectoralwin ~ treat_vign, data=f[as.numeric(f$treat_vign)==1 | as.numeric(f$treat_vign)==i,])$p.value
}
p.adjust(pvallist,method="BH")



# ================================================================
# Appendix A6: Beliefs about Germany by Leavers/Remainers
# ================================================================

f$BrexitBin <- ifelse(f$Brexit==1,1,ifelse(f$Brexit==2,2,NA))
subdf <- f[is.na(f$BrexitBin)==FALSE,]
subdf$BrexitBin <- factor(subdf$BrexitBin, levels=c(1,2), labels=c("Leavers", "Remainers"))

# Beliefs for Leavers
h1 <- group_by(subdf[subdf$BrexitBin=="Leavers",], treat_vign) %>%
  summarise(
    count = n(),
    mean = mean(vign_beliefGerwins, na.rm = TRUE),
    sd = sd(vign_beliefGerwins, na.rm = TRUE)
  )


# Beliefs for Remainers
h2 <- group_by(subdf[subdf$BrexitBin=="Remainers",], treat_vign) %>%
  summarise(
    count = n(),
    mean = mean(vign_beliefGerwins, na.rm = TRUE),
    sd = sd(vign_beliefGerwins, na.rm = TRUE)
  )


# Compare Leavers vs Remainers
# Multiple comparison p-value correction
pvallist <- NA
for (i in 1:5) {
  pvallist[i] <- t.test(vign_beliefGerwins ~ BrexitBin, data=subdf[as.numeric(subdf$treat_vign)==i,])$p.value
}
p.adjust(pvallist,method="BH")


# ================================================================
# Figure A3: Leavers vs Remainers
# ================================================================


# Restructure data for better plotting options
subdf.plot <- cbind(subdf[,c("treat_vign","BrexitBin","vign_beliefGerwins")])
colnames(subdf.plot)[2:3] <- c("variable", "outcome")


pdf(file="./belief_Ger_win_Brexit.pdf", width=8, height=5)

cols <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")

# Construct data set for means
meandf <- data.frame(variable=c(levels(subdf.plot$variable)), mean=c(h1$mean[1],h2$mean[1]))

ggerrorplot(subdf.plot, x="treat_vign", y="outcome", color="variable",
            desc_stat="mean_ci", 
            error.plot="errorbar",
            add="mean",
            legend="none",
            ylab="DV (ordered)", xlab = "Experimental conditions",
            palette=c(cols[2],cols[3]),
            ylim=c(1.5,3.5)
) +
  facet_grid(~variable, scales="free") +
  scale_x_discrete(breaks=c("C", "Tr_DomWin", "Tr_DomLose", "Tr_ForWin", "Tr_ForLose"), labels=c("Control", "T1: UK \n sector wins", "T2: UK \n sector loses", "T3: German \n sector wins", "T4: German \n sector loses")) +
  geom_hline(data=meandf, aes(yintercept=mean), linetype=c("dashed","dashed"), size=.5,col=c(cols[2],cols[3])) +
  theme_bw() +
  theme(text = element_text(size=11),
        axis.text.x = element_text(size=9),
        legend.position="none"
  )
dev.off()

# ================================================================
# Appendix A7: Heterogeneous Treatment Effects
# ================================================================

# ================================================================
# Appendix A7.1: Sector Employment
# ================================================================


f$secemploy <- ifelse((f$empl_sector==3 | f$empl_sector==4 | f$empl_sector==7 | f$empl_sector==9),1,ifelse((f$empl_sector==1 | f$empl_sector==2 | f$empl_sector==5 | f$empl_sector==6 | f$empl_sector==8 | f$empl_sector==10),2,NA))
 

# Restructure data for better plotting options
subdf <- f[is.na(f$secemploy)==FALSE,]
subdf$secemploy <- factor(subdf$secemploy, levels=c(1,2), labels=c("Highly carbon intensive sector", "Low carbon intensive sector"))


# Beliefs for carbon intensive sectors
h1 <- group_by(subdf[subdf$secemploy=="Highly carbon intensive sector",], treat_vign) %>%
  summarise(
    count = n(),
    mean = mean(vign_beliefUKwins, na.rm = TRUE),
    sd = sd(vign_beliefUKwins, na.rm = TRUE)
  )

# ANOVA: Highly carbon intensive sectors
res.aov <- aov(vign_beliefUKwins ~ treat_vign, data = subdf[subdf$secemploy=="Highly carbon intensive sector",])
summary(res.aov)
TukeyHSD(res.aov)

# Beliefs for low carbon intensive sectors
h2 <- group_by(subdf[subdf$secemploy=="Low carbon intensive sector",], treat_vign) %>%
  summarise(
    count = n(),
    mean = mean(vign_beliefUKwins, na.rm = TRUE),
    sd = sd(vign_beliefUKwins, na.rm = TRUE)
  )

# ANOVA: Low carbon intensive sectors
res.aov <- aov(vign_beliefUKwins ~ treat_vign, data = subdf[subdf$secemploy=="Low carbon intensive sector",])
summary(res.aov)
TukeyHSD(res.aov)


# Compare subgroups
# Multiple comparison p-value correction
pvallist <- NA
for (i in 1:5) {
  pvallist[i] <- t.test(vign_beliefUKwins ~ secemploy, data=subdf[as.numeric(subdf$treat_vign)==i,])$p.value
}
p.adjust(pvallist,method="BH")


# ================================================================
# Figure A4: Highly carbon intensive sectors vs low carbon intensive sectors
# ================================================================

# Restructure data for better plotting options
subdf.plot <- cbind(subdf[,c("treat_vign","secemploy","vign_beliefUKwins")])
colnames(subdf.plot)[2:3] <- c("variable", "outcome")


pdf(file="./belief_UK_win_secemploy.pdf", width=8, height=5)

cols <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")

# Construct data set for means
meandf <- data.frame(variable=c(levels(subdf.plot$variable)), mean=c(h1$mean[1],h2$mean[1]))

ggerrorplot(subdf.plot, x="treat_vign", y="outcome", color="variable",
            desc_stat="mean_ci", 
            error.plot="errorbar",
            add="mean",
            legend="none",
            ylab="DV (ordered)", xlab = "Experimental conditions",
            palette=c(cols[2],cols[3]),
            ylim=c(1.5,3.5)
) +
  facet_grid(~variable, scales="free") +
  scale_x_discrete(breaks=c("C", "Tr_DomWin", "Tr_DomLose", "Tr_ForWin", "Tr_ForLose"), labels=c("Control", "T1: UK \n sector wins", "T2: UK \n sector loses", "T3: German \n sector wins", "T4: German \n sector loses")) +
  geom_hline(data=meandf, aes(yintercept=mean), linetype=c("dashed","dashed"), size=.5,col=c(cols[2],cols[3])) +
  theme_bw() +
  theme(text = element_text(size=11),
        axis.text.x = element_text(size=9),
        legend.position="none"
  )
dev.off()



# ================================================================
# Appendix A7.2 Extractive Industries
# ================================================================

# (1) COAL MINING REGIONS

f$geo <- ifelse((f$residence==2 | f$residence==3),1,ifelse(f$gor_england<6,1,2))
subdf <- f[is.na(f$geo)==FALSE,]
subdf$geo <- factor(subdf$geo, levels=c(1,2), labels=c("Coal mining region", "No coal mining"))


# Beliefs for coal regions
h1 <- group_by(subdf[subdf$geo=="Coal mining region",], treat_vign) %>%
  summarise(
    count = n(),
    mean = mean(vign_beliefUKwins, na.rm = TRUE),
    sd = sd(vign_beliefUKwins, na.rm = TRUE)
  )

# ANOVA: coal regions
res.aov <- aov(vign_beliefUKwins ~ treat_vign, data = subdf[subdf$geo=="Coal mining region",])
summary(res.aov)
TukeyHSD(res.aov)

# Beliefs for non-coal regions
h2 <- group_by(subdf[subdf$geo=="No coal mining",], treat_vign) %>%
  summarise(
    count = n(),
    mean = mean(vign_beliefUKwins, na.rm = TRUE),
    sd = sd(vign_beliefUKwins, na.rm = TRUE)
  )

# ANOVA: Non-coal regions
res.aov <- aov(vign_beliefUKwins ~ treat_vign, data = subdf[subddf$geo=="No coal mining",])
summary(res.aov)
TukeyHSD(res.aov)

# Compare subgroups
# Multiple comparison p-value correction
pvallist <- NA
for (i in 1:5) {
  pvallist[i] <- t.test(vign_beliefUKwins ~ geo, data=subdf[as.numeric(subdf$treat_vign)==i,])$p.value
}
p.adjust(pvallist,method="BH")


# ================================================================
# Figure A5: Coal mining region vs no coal mining region
# ================================================================


# Restructure data for better plotting options
subdf.plot <- cbind(subdf[,c("treat_vign","geo","vign_beliefUKwins")])
colnames(subdf.plot)[2:3] <- c("variable", "outcome")


pdf(file="./belief_UK_win_coal.pdf", width=8, height=5)

cols <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")

# Construct data set for means
meandf <- data.frame(variable=c(levels(subdf.plot$variable)), mean=c(h1$mean[1],h2$mean[1]))

ggerrorplot(subdf.plot, x="treat_vign", y="outcome", color="variable",
            desc_stat="mean_ci", 
            error.plot="errorbar",
            add="mean",
            legend="none",
            ylab="DV (ordered)", xlab = "Experimental conditions",
            palette=c(cols[2],cols[3]),
            ylim=c(1.5,3.5)
) +
  facet_grid(~variable, scales="free") +
  scale_x_discrete(breaks=c("C", "Tr_DomWin", "Tr_DomLose", "Tr_ForWin", "Tr_ForLose"), labels=c("Control", "T1: UK \n sector wins", "T2: UK \n sector loses", "T3: German \n sector wins", "T4: German \n sector loses")) +
  geom_hline(data=meandf, aes(yintercept=mean), linetype=c("dashed","dashed"), size=.5,col=c(cols[2],cols[3])) +
  theme_bw() +
  theme(text = element_text(size=11),
        axis.text.x = element_text(size=9),
        legend.position="none"
  )
dev.off()



# (2) OIL REGIONS

f$geo <- ifelse(f$residence==3,1,ifelse((f$gor_england==1 | f$gor_england==2 | f$gor_england==4 | f$gor_england==6),1,2))
subdf <- f[is.na(f$geo)==FALSE,]
subdf$geo <- factor(subdf$geo, levels=c(1,2), labels=c("Oil/gas producing region", "No oil/gas production"))


# Beliefs for oil/gas regions
h1 <- group_by(subdf[subdf$geo=="Oil/gas producing region",], treat_vign) %>%
  summarise(
    count = n(),
    mean = mean(vign_beliefUKwins, na.rm = TRUE),
    sd = sd(vign_beliefUKwins, na.rm = TRUE)
  )

# ANOVA: Oil/gas regions
res.aov <- aov(vign_beliefUKwins ~ treat_vign, data = subdf[subdf$geo=="Oil/gas producing region",])
summary(res.aov)
TukeyHSD(res.aov)

# Beliefs for no oil/gas region
h2 <- group_by(subdf[subdf$geo=="No oil/gas production",], treat_vign) %>%
  summarise(
    count = n(),
    mean = mean(vign_beliefUKwins, na.rm = TRUE),
    sd = sd(vign_beliefUKwins, na.rm = TRUE)
  )

# ANOVA: No oil/gas regions
res.aov <- aov(vign_beliefUKwins ~ treat_vign, data = subdf[subddf$geo=="No oil/gas production",])
summary(res.aov)
TukeyHSD(res.aov)

# Compare subgroups
# Multiple comparison p-value correction
pvallist <- NA
for (i in 1:5) {
  pvallist[i] <- t.test(vign_beliefUKwins ~ geo, data=subdf[as.numeric(subdf$treat_vign)==i,])$p.value
}
p.adjust(pvallist,method="BH")


# ================================================================
# Figure A6: Oil/gas producing region vs no oil gas production
# ================================================================


# Restructure data for better plotting options
subdf.plot <- cbind(subdf[,c("treat_vign","geo","vign_beliefUKwins")])
colnames(subdf.plot)[2:3] <- c("variable", "outcome")


pdf(file="./belief_UK_win_oil.pdf", width=8, height=5)

cols <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")

# Construct data set for means
meandf <- data.frame(variable=c(levels(subdf.plot$variable)), mean=c(h1$mean[1],h2$mean[1]))

ggerrorplot(subdf.plot, x="treat_vign", y="outcome", color="variable",
            desc_stat="mean_ci", 
            error.plot="errorbar",
            add="mean",
            legend="none",
            ylab="DV (ordered)", xlab = "Experimental conditions",
            palette=c(cols[2],cols[3]),
            ylim=c(1.5,3.5)
) +
  facet_grid(~variable, scales="free") +
  scale_x_discrete(breaks=c("C", "Tr_DomWin", "Tr_DomLose", "Tr_ForWin", "Tr_ForLose"), labels=c("Control", "T1: UK \n sector wins", "T2: UK \n sector loses", "T3: German \n sector wins", "T4: German \n sector loses")) +
  geom_hline(data=meandf, aes(yintercept=mean), linetype=c("dashed","dashed"), size=.5,col=c(cols[2],cols[3])) +
  theme_bw() +
  theme(text = element_text(size=11),
        axis.text.x = element_text(size=9),
        legend.position="none"
  )
dev.off()



# ================================================================
# Appendix A7.3 Geography
# ================================================================


# (1) ENGLAND VERSUS SCOTLAND

f$geo <- ifelse(f$residence==1,1,ifelse(f$residence==3,2,NA))
subdf <- f[is.na(f$geo)==FALSE,]
subdf$geo <- factor(subdf$geo, levels=c(1,2), labels=c("England", "Scotland"))


# Beliefs for England
h1 <- group_by(subdf[subdf$geo=="England",], treat_vign) %>%
  summarise(
    count = n(),
    mean = mean(vign_beliefUKwins, na.rm = TRUE),
    sd = sd(vign_beliefUKwins, na.rm = TRUE)
  )

# ANOVA: England
res.aov <- aov(vign_beliefUKwins ~ treat_vign, data = subdf[subdf$geo=="England",])
summary(res.aov)
TukeyHSD(res.aov)

# Beliefs for Scotland
h2 <- group_by(subdf[subdf$geo=="Scotland",], treat_vign) %>%
  summarise(
    count = n(),
    mean = mean(vign_beliefUKwins, na.rm = TRUE),
    sd = sd(vign_beliefUKwins, na.rm = TRUE)
  )

# ANOVA: Scotland
res.aov <- aov(vign_beliefUKwins ~ treat_vign, data = subdf[subdf$geo=="Scotland",])
summary(res.aov)
TukeyHSD(res.aov)

# Compare subgroups
# Multiple comparison p-value correction
pvallist <- NA
for (i in 1:5) {
  pvallist[i] <- t.test(vign_beliefUKwins ~ geo, data=subdf[as.numeric(subdf$treat_vign)==i,])$p.value
}
p.adjust(pvallist,method="BH")


# ================================================================
# Figure A7: England vs Scotland
# ================================================================


# Restructure data for better plotting options
subdf.plot <- cbind(subdf[,c("treat_vign","geo","vign_beliefUKwins")])
colnames(subdf.plot)[2:3] <- c("variable", "outcome")


pdf(file="./belief_UK_win_EngScot.pdf", width=8, height=5)

cols <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")

# Construct data set for means
meandf <- data.frame(variable=c(levels(subdf.plot$variable)), mean=c(h1$mean[1],h2$mean[1]))

ggerrorplot(subdf.plot, x="treat_vign", y="outcome", color="variable",
            desc_stat="mean_ci", 
            error.plot="errorbar",
            add="mean",
            legend="none",
            ylab="DV (ordered)", xlab = "Experimental conditions",
            palette=c(cols[2],cols[3]),
            ylim=c(1.5,3.5)
) +
  facet_grid(~variable, scales="free") +
  scale_x_discrete(breaks=c("C", "Tr_DomWin", "Tr_DomLose", "Tr_ForWin", "Tr_ForLose"), labels=c("Control", "T1: UK \n sector wins", "T2: UK \n sector loses", "T3: German \n sector wins", "T4: German \n sector loses")) +
  geom_hline(data=meandf, aes(yintercept=mean), linetype=c("dashed","dashed"), size=.5,col=c(cols[2],cols[3])) +
  theme_bw() +
  theme(text = element_text(size=11),
        axis.text.x = element_text(size=9),
        legend.position="none"
  )
dev.off()



# (2) LONDON VERSUS NON-LONDON

f$geo <- ifelse(f$gor_england==7,1,ifelse(f$gor_england!=7,2,NA))
subdf <- f[is.na(f$geo)==FALSE,]
subdf$geo <- factor(subdf$geo, levels=c(1,2), labels=c("London", "Outside London"))


# Beliefs for London
h1 <- group_by(subdf[subdf$geo=="London",], treat_vign) %>%
  summarise(
    count = n(),
    mean = mean(vign_beliefUKwins, na.rm = TRUE),
    sd = sd(vign_beliefUKwins, na.rm = TRUE)
  )

# ANOVA: London
res.aov <- aov(vign_beliefUKwins ~ treat_vign, data = subdf[subdf$geo=="London",])
summary(res.aov)
TukeyHSD(res.aov)

# Beliefs for Outside London
h2 <- group_by(subdf[subdf$geo=="Outside London",], treat_vign) %>%
  summarise(
    count = n(),
    mean = mean(vign_beliefUKwins, na.rm = TRUE),
    sd = sd(vign_beliefUKwins, na.rm = TRUE)
  )


# ANOVA: Outside London
res.aov <- aov(vign_beliefUKwins ~ treat_vign, data = subdf[subdf$geo=="Outside London",])
summary(res.aov)
TukeyHSD(res.aov)

# Compare subgroups
# Multiple comparison p-value correction
pvallist <- NA
for (i in 1:5) {
  pvallist[i] <- t.test(vign_beliefUKwins ~ geo, data=subdf[as.numeric(subdf$treat_vign)==i,])$p.value
}
p.adjust(pvallist,method="BH")


# ================================================================
# Figure A8: London vs outside London
# ================================================================


# Restructure data for better plotting options
subdf.plot <- cbind(subdf[,c("treat_vign","geo","vign_beliefUKwins")])
colnames(subdf.plot)[2:3] <- c("variable", "outcome")


pdf(file="./belief_UK_win_London.pdf", width=8, height=5)

cols <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")

# Construct data set for means
meandf <- data.frame(variable=c(levels(subdf.plot$variable)), mean=c(h1$mean[1],h2$mean[1]))

ggerrorplot(subdf.plot, x="treat_vign", y="outcome", color="variable",
            desc_stat="mean_ci", 
            error.plot="errorbar",
            add="mean",
            legend="none",
            ylab="DV (ordered)", xlab = "Experimental conditions",
            palette=c(cols[2],cols[3]),
            ylim=c(1.5,3.5)
) +
  facet_grid(~variable, scales="free") +
  scale_x_discrete(breaks=c("C", "Tr_DomWin", "Tr_DomLose", "Tr_ForWin", "Tr_ForLose"), labels=c("Control", "T1: UK \n sector wins", "T2: UK \n sector loses", "T3: German \n sector wins", "T4: German \n sector loses")) +
  geom_hline(data=meandf, aes(yintercept=mean), linetype=c("dashed","dashed"), size=.5,col=c(cols[2],cols[3])) +
  theme_bw() +
  theme(text = element_text(size=11),
        axis.text.x = element_text(size=9),
        legend.position="none"
  )
dev.off()


# (3) NORTH SOUTH DIVIDE

f$geo <- ifelse(f$gor_england>=6,1,ifelse(f$gor_england<6,2,NA))
subdf <- f[is.na(f$geo)==FALSE,]
subdf$geo <- factor(subdf$geo, levels=c(1,2), labels=c("South", "North"))


# Beliefs for South
h1 <- group_by(subdf[subdf$geo=="South",], treat_vign) %>%
  summarise(
    count = n(),
    mean = mean(vign_beliefUKwins, na.rm = TRUE),
    sd = sd(vign_beliefUKwins, na.rm = TRUE)
  )

# ANOVA: South
res.aov <- aov(vign_beliefUKwins ~ treat_vign, data = subdf[subdf$geo=="South",])
summary(res.aov)
TukeyHSD(res.aov)

# Beliefs for North
h2 <- group_by(subdf[subdf$geo=="North",], treat_vign) %>%
  summarise(
    count = n(),
    mean = mean(vign_beliefUKwins, na.rm = TRUE),
    sd = sd(vign_beliefUKwins, na.rm = TRUE)
  )

# ANOVA: Outside North
res.aov <- aov(vign_beliefUKwins ~ treat_vign, data = subdf[subddf$geo=="North",])
summary(res.aov)
TukeyHSD(res.aov)

# Compare subgroups
# Multiple comparison p-value correction
pvallist <- NA
for (i in 1:5) {
  pvallist[i] <- t.test(vign_beliefUKwins ~ geo, data=subdf[as.numeric(subdf$treat_vign)==i,])$p.value
}
p.adjust(pvallist,method="BH")


# ================================================================
# Figure A9: North vs South
# ================================================================


# Restructure data for better plotting options
subdf.plot <- cbind(subdf[,c("treat_vign","geo","vign_beliefUKwins")])
colnames(subdf.plot)[2:3] <- c("variable", "outcome")


pdf(file="./belief_UK_win_NoSo.pdf", width=8, height=5)

cols <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")

# Construct data set for means
meandf <- data.frame(variable=c(levels(subdf.plot$variable)), mean=c(h1$mean[1],h2$mean[1]))

ggerrorplot(subdf.plot, x="treat_vign", y="outcome", color="variable",
            desc_stat="mean_ci", 
            error.plot="errorbar",
            add="mean",
            legend="none",
            ylab="DV (ordered)", xlab = "Experimental conditions",
            palette=c(cols[2],cols[3]),
            ylim=c(1.5,3.5)
) +
  facet_grid(~variable, scales="free") +
  scale_x_discrete(breaks=c("C", "Tr_DomWin", "Tr_DomLose", "Tr_ForWin", "Tr_ForLose"), labels=c("Control", "T1: UK \n sector wins", "T2: UK \n sector loses", "T3: German \n sector wins", "T4: German \n sector loses")) +
  geom_hline(data=meandf, aes(yintercept=mean), linetype=c("dashed","dashed"), size=.5,col=c(cols[2],cols[3])) +
  theme_bw() +
  theme(text = element_text(size=11),
        axis.text.x = element_text(size=9),
        legend.position="none"
  )
dev.off()



# ================================================================
# Appendix A7.4: Income
# ================================================================

f$inc <- ifelse(f$income>3 & f$income<9,1,ifelse(f$sectorperformance<4,2,NA))

# Restructure data for better plotting options
subdf <- f[is.na(f$inc)==FALSE,]
subdf$inc <- factor(subdf$inc, levels=c(1,2), labels=c("High income", "Low income"))

# Beliefs for high income
h1 <- group_by(subdf[subdf$inc=="High income",], treat_vign) %>%
  summarise(
    count = n(),
    mean = mean(vign_beliefUKwins, na.rm = TRUE),
    sd = sd(vign_beliefUKwins, na.rm = TRUE)
  )

# Beliefs for low income
h2 <- group_by(subdf[subdf$inc=="Low income",], treat_vign) %>%
  summarise(
    count = n(),
    mean = mean(vign_beliefUKwins, na.rm = TRUE),
    sd = sd(vign_beliefUKwins, na.rm = TRUE)
  )


# Compare subgroups
# Multiple comparison p-value correction
pvallist <- NA
for (i in 1:5) {
  pvallist[i] <- t.test(vign_beliefUKwins ~ inc, data=subdf[as.numeric(subdf$treat_vign)==i,])$p.value
}
p.adjust(pvallist,method="BH")


# ================================================================
# Figure A10: Income
# ================================================================

# Restructure data for better plotting options
subdf.plot <- cbind(subdf[,c("treat_vign","inc","vign_beliefUKwins")])
colnames(subdf.plot)[2:3] <- c("variable", "outcome")


pdf(file="./belief_UK_win_income.pdf", width=8, height=5)

cols <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")

# Construct data set for means
meandf <- data.frame(variable=c(levels(subdf.plot$variable)), mean=c(h1$mean[1],h2$mean[1]))

ggerrorplot(subdf.plot, x="treat_vign", y="outcome", color="variable",
            desc_stat="mean_ci", 
            error.plot="errorbar",
            add="mean",
            legend="none",
            ylab="DV (ordered)", xlab = "Experimental conditions",
            palette=c(cols[2],cols[3]),
            ylim=c(1.5,3.5)
) +
  facet_grid(~variable, scales="free") +
  scale_x_discrete(breaks=c("C", "Tr_DomWin", "Tr_DomLose", "Tr_ForWin", "Tr_ForLose"), labels=c("Control", "T1: UK \n sector wins", "T2: UK \n sector loses", "T3: German \n sector wins", "T4: German \n sector loses")) +
  geom_hline(data=meandf, aes(yintercept=mean), linetype=c("dashed","dashed"), size=.5,col=c(cols[2],cols[3])) +
  theme_bw() +
  theme(text = element_text(size=11),
        axis.text.x = element_text(size=9),
        legend.position="none"
  )
dev.off()



# ================================================================
# Appendix A7.5: Party Identification
# ================================================================

# (1) PARTY IDENTIFICATION

f$pol_id <- ifelse(f$political_id<=4,1,ifelse(f$political_id>4,2,NA))

# Restructure data for better plotting options
subdf <- f[is.na(f$pol_id)==FALSE,]
subdf$pol_id <- factor(subdf$pol_id, levels=c(1,2), labels=c("Left/Center", "Right"))


# Beliefs for Leavers
h1 <- group_by(subdf[subdf$pol_id=="Left/Center",], treat_vign) %>%
  summarise(
    count = n(),
    mean = mean(vign_beliefUKwins, na.rm = TRUE),
    sd = sd(vign_beliefUKwins, na.rm = TRUE)
  )

# Beliefs for Remainers
h2 <- group_by(subdf[subdf$pol_id=="Right",], treat_vign) %>%
  summarise(
    count = n(),
    mean = mean(vign_beliefUKwins, na.rm = TRUE),
    sd = sd(vign_beliefUKwins, na.rm = TRUE)
  )


# Compare subgroups
# Multiple comparison p-value correction
pvallist <- NA
for (i in 1:5) {
  pvallist[i] <- t.test(vign_beliefUKwins ~ pol_id, data=subdf[as.numeric(subdf$treat_vign)==i,])$p.value
}
p.adjust(pvallist,method="BH")


# ================================================================
# Figure A11: Left/centre vs right
# ================================================================


# Restructure data for better plotting options
subdf.plot <- cbind(subdf[,c("treat_vign","pol_id","vign_beliefUKwins")])
colnames(subdf.plot)[2:3] <- c("variable", "outcome")


pdf(file="./belief_UK_win_polid.pdf", width=8, height=5)

cols <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")

# Construct data set for means
meandf <- data.frame(variable=c(levels(subdf.plot$variable)), mean=c(h1$mean[1],h2$mean[1]))

ggerrorplot(subdf.plot, x="treat_vign", y="outcome", color="variable",
            desc_stat="mean_ci", 
            error.plot="errorbar",
            add="mean",
            legend="none",
            ylab="DV (ordered)", xlab = "Experimental conditions",
            palette=c(cols[2],cols[3]),
            ylim=c(1.5,3.5)
) +
  facet_grid(~variable, scales="free") +
  scale_x_discrete(breaks=c("C", "Tr_DomWin", "Tr_DomLose", "Tr_ForWin", "Tr_ForLose"), labels=c("Control", "T1: UK \n sector wins", "T2: UK \n sector loses", "T3: German \n sector wins", "T4: German \n sector loses")) +
  geom_hline(data=meandf, aes(yintercept=mean), linetype=c("dashed","dashed"), size=.5,col=c(cols[2],cols[3])) +
  theme_bw() +
  theme(text = element_text(size=11),
        axis.text.x = element_text(size=9),
        legend.position="none"
  )
dev.off()


# (2) VOTE CHOICE

f$elecparty <- ifelse(f$parties==1,1,ifelse(f$parties==2,2,NA))

# Restructure data for better plotting options
df.plot <- f[is.na(f$elecparty)==FALSE,]
df.plot$elecparty <- factor(df.plot$elecparty, levels=c(1,2), labels=c("Tories", "Labour"))


# Beliefs for Tory voters
group_by(df.plot[df.plot$elecparty=="Tories",], treat_vign) %>%
  summarise(
    count = n(),
    mean = mean(vign_beliefUKwins, na.rm = TRUE),
    sd = sd(vign_beliefUKwins, na.rm = TRUE)
  )

# Beliefs for Labour voters
group_by(df.plot[df.plot$elecparty=="Labour",], treat_vign) %>%
  summarise(
    count = n(),
    mean = mean(vign_beliefUKwins, na.rm = TRUE),
    sd = sd(vign_beliefUKwins, na.rm = TRUE)
  )



# Compare respondents across groups
# Multiple comparison p-value correction
pvallist <- NA
for (i in 1:5) {
  pvallist[i] <- t.test(vign_beliefUKwins ~ elecparty, data=df.plot[as.numeric(df.plot$treat_vign)==i ,])$p.value
}
p.adjust(pvallist,method="BH")



# ================================================================
# Appendix A7.6: Sector Attitudes
# ================================================================

# (1) SECTOR PERFORMANCE

f$sec <- ifelse(f$sectorperformance>2,1,ifelse(f$sectorperformance<3,2,NA))

# Restructure data for better plotting options
subdf <- f[is.na(f$sec)==FALSE,]
subdf$sec <- factor(subdf$sec, levels=c(1,2), labels=c("Sectors do well", "Sectors do badly"))


# Beliefs for sectors do well
h1 <- group_by(subdf[subdf$sec=="Sectors do well",], treat_vign) %>%
  summarise(
    count = n(),
    mean = mean(vign_beliefUKwins, na.rm = TRUE),
    sd = sd(vign_beliefUKwins, na.rm = TRUE)
  )

# Beliefs for Remainers
h2 <- group_by(subdf[subdf$sec=="Sectors do badly",], treat_vign) %>%
  summarise(
    count = n(),
    mean = mean(vign_beliefUKwins, na.rm = TRUE),
    sd = sd(vign_beliefUKwins, na.rm = TRUE)
  )


# Compare subgroups
# Multiple comparison p-value correction
pvallist <- NA
for (i in 1:5) {
  pvallist[i] <- t.test(vign_beliefUKwins ~ sec, data=subdf[as.numeric(subdf$treat_vign)==i,])$p.value
}
p.adjust(pvallist,method="BH")


# ================================================================
# Figure A12: Sectors do well/badly
# ================================================================

# Restructure data for better plotting options
subdf.plot <- cbind(subdf[,c("treat_vign","sec","vign_beliefUKwins")])
colnames(subdf.plot)[2:3] <- c("variable", "outcome")


pdf(file="./belief_UK_win_sector.pdf", width=8, height=5)

cols <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")

# Construct data set for means
meandf <- data.frame(variable=c(levels(subdf.plot$variable)), mean=c(h1$mean[1],h2$mean[1]))

ggerrorplot(subdf.plot, x="treat_vign", y="outcome", color="variable",
            desc_stat="mean_ci", 
            error.plot="errorbar",
            add="mean",
            legend="none",
            ylab="DV (ordered)", xlab = "Experimental conditions",
            palette=c(cols[2],cols[3]),
            ylim=c(1.5,3.5)
) +
  facet_grid(~variable, scales="free") +
  scale_x_discrete(breaks=c("C", "Tr_DomWin", "Tr_DomLose", "Tr_ForWin", "Tr_ForLose"), labels=c("Control", "T1: UK \n sector wins", "T2: UK \n sector loses", "T3: German \n sector wins", "T4: German \n sector loses")) +
  geom_hline(data=meandf, aes(yintercept=mean), linetype=c("dashed","dashed"), size=.5,col=c(cols[2],cols[3])) +
  theme_bw() +
  theme(text = element_text(size=11),
        axis.text.x = element_text(size=9),
        legend.position="none"
  )
dev.off()


# (2) PUBLIC TRANSPORT USE

f$trans <- ifelse(f$transportsectorfeel==1,1,ifelse(f$transportsectorfeel==4,2,NA))

# Restructure data for better plotting options
df.plot <- f[is.na(f$trans)==FALSE,]
df.plot$trans <- factor(df.plot$trans, levels=c(1,2), labels=c("Public transport", "No public transport"))


# Beliefs for: Sectors do well
group_by(df.plot[df.plot$trans=="Public transport",], treat_vign) %>%
  summarise(
    count = n(),
    mean = mean(vign_beliefUKwins, na.rm = TRUE),
    sd = sd(vign_beliefUKwins, na.rm = TRUE)
  )

# Beliefs for: Sectors do badly
group_by(df.plot[df.plot$trans=="No public transport",], treat_vign) %>%
  summarise(
    count = n(),
    mean = mean(vign_beliefUKwins, na.rm = TRUE),
    sd = sd(vign_beliefUKwins, na.rm = TRUE)
  )

# Compare respondents across groups
# Multiple comparison p-value correction
pvallist <- NA
for (i in 1:5) {
  pvallist[i] <- t.test(vign_beliefUKwins ~ trans, data=df.plot[as.numeric(df.plot$treat_vign)==i ,])$p.value
}
p.adjust(pvallist,method="BH")

# (3) UNDERSTAND GAS/ELECTRICITY BILL

f$energy <- ifelse(f$energysectorfeel<3,1,ifelse(f$energysectorfeel>2,2,NA))

# Restructure data for better plotting options
df.plot <- f[is.na(f$energy)==FALSE,]
df.plot$energy <- factor(df.plot$energy, levels=c(1,2), labels=c("High knowledge", "Little knowledge"))


# Beliefs for: Understand electricity bill well
group_by(df.plot[df.plot$energy=="High knowledge",], treat_vign) %>%
  summarise(
    count = n(),
    mean = mean(vign_beliefUKwins, na.rm = TRUE),
    sd = sd(vign_beliefUKwins, na.rm = TRUE)
  )

# Beliefs for: Do not understand the electricity bill well
group_by(df.plot[df.plot$energy=="Little knowledge",], treat_vign) %>%
  summarise(
    count = n(),
    mean = mean(vign_beliefUKwins, na.rm = TRUE),
    sd = sd(vign_beliefUKwins, na.rm = TRUE)
  )

# Compare respondents across groups
# Multiple comparison p-value correction
pvallist <- NA
for (i in 1:5) {
  pvallist[i] <- t.test(vign_beliefUKwins ~ energy, data=df.plot[as.numeric(df.plot$treat_vign)==i ,])$p.value
}
p.adjust(pvallist,method="BH")



# ================================================================
# Appendix A7.7: Climate Attitude
# ================================================================

# (1) CLIMATE CHANGE IS A SERIOUS PROBLEM

f$serious <- ifelse(f$cc_seriousproblem<2,1,ifelse(f$cc_seriousproblem>1,2,NA))

# Restructure data for better plotting options
subdf <- f[is.na(f$serious)==FALSE,]
subdf$serious <- factor(subdf$serious, levels=c(1,2), labels=c("Very serious problem", "Not so serious problem"))


# Beliefs for CC is a serious problem
h1 <- group_by(subdf[subdf$serious=="Very serious problem",], treat_vign) %>%
  summarise(
    count = n(),
    mean = mean(vign_beliefUKwins, na.rm = TRUE),
    sd = sd(vign_beliefUKwins, na.rm = TRUE)
  )

# Beliefs for Remainers
h2 <- group_by(subdf[subdf$serious=="Not so serious problem",], treat_vign) %>%
  summarise(
    count = n(),
    mean = mean(vign_beliefUKwins, na.rm = TRUE),
    sd = sd(vign_beliefUKwins, na.rm = TRUE)
  )


# Compare subgroups
# Multiple comparison p-value correction
pvallist <- NA
for (i in 1:5) {
  pvallist[i] <- t.test(vign_beliefUKwins ~ serious, data=subdf[as.numeric(subdf$treat_vign)==i,])$p.value
}
p.adjust(pvallist,method="BH")


# ================================================================
# Figure A13: Very serious problem vs not so serious problem
# ================================================================


# Restructure data for better plotting options
subdf.plot <- cbind(subdf[,c("treat_vign","serious","vign_beliefUKwins")])
colnames(subdf.plot)[2:3] <- c("variable", "outcome")


pdf(file="./belief_UK_win_ccproblem.pdf", width=8, height=5)

cols <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")

# Construct data set for means
meandf <- data.frame(variable=c(levels(subdf.plot$variable)), mean=c(h1$mean[1],h2$mean[1]))

ggerrorplot(subdf.plot, x="treat_vign", y="outcome", color="variable",
            desc_stat="mean_ci", 
            error.plot="errorbar",
            add="mean",
            legend="none",
            ylab="DV (ordered)", xlab = "Experimental conditions",
            palette=c(cols[2],cols[3]),
            ylim=c(1.5,3.5)
) +
  facet_grid(~variable, scales="free") +
  scale_x_discrete(breaks=c("C", "Tr_DomWin", "Tr_DomLose", "Tr_ForWin", "Tr_ForLose"), labels=c("Control", "T1: UK \n sector wins", "T2: UK \n sector loses", "T3: German \n sector wins", "T4: German \n sector loses")) +
  geom_hline(data=meandf, aes(yintercept=mean), linetype=c("dashed","dashed"), size=.5,col=c(cols[2],cols[3])) +
  theme_bw() +
  theme(text = element_text(size=11),
        axis.text.x = element_text(size=9),
        legend.position="none"
  )
dev.off()


# (2) CONCERN ABOUT CLIMATE CHANGE

f$concern <- ifelse(f$cc_concern>2,1,ifelse(f$cc_concern<3,2,NA))

# Restructure data for better plotting options
subdf <- f[is.na(f$concern)==FALSE,]
subdf$concern <- factor(subdf$concern, levels=c(1,2), labels=c("Very concerned", "Not so concerned"))


# Beliefs for those concerned about CC
h1 <- group_by(subdf[subdf$concern=="Very concerned",], treat_vign) %>%
  summarise(
    count = n(),
    mean = mean(vign_beliefUKwins, na.rm = TRUE),
    sd = sd(vign_beliefUKwins, na.rm = TRUE)
  )

# Beliefs for Remainers
h2 <- group_by(subdf[subdf$concern=="Not so concerned",], treat_vign) %>%
  summarise(
    count = n(),
    mean = mean(vign_beliefUKwins, na.rm = TRUE),
    sd = sd(vign_beliefUKwins, na.rm = TRUE)
  )


# Compare Leavers vs Remainers
# Multiple comparison p-value correction
pvallist <- NA
for (i in 1:5) {
  pvallist[i] <- t.test(vign_beliefUKwins ~ concern, data=subdf[as.numeric(subdf$treat_vign)==i,])$p.value
}
p.adjust(pvallist,method="BH")


# ================================================================
# Figure A14: Very concerned vs not concerned
# ================================================================


# Restructure data for better plotting options
subdf.plot <- cbind(subdf[,c("treat_vign","concern","vign_beliefUKwins")])
colnames(subdf.plot)[2:3] <- c("variable", "outcome")


pdf(file="./belief_UK_win_ccconcern.pdf", width=8, height=5)

cols <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")

# Construct data set for means
meandf <- data.frame(variable=c(levels(subdf.plot$variable)), mean=c(h1$mean[1],h2$mean[1]))

ggerrorplot(subdf.plot, x="treat_vign", y="outcome", color="variable",
            desc_stat="mean_ci", 
            error.plot="errorbar",
            add="mean",
            legend="none",
            ylab="DV (ordered)", xlab = "Experimental conditions",
            palette=c(cols[2],cols[3]),
            ylim=c(1.5,3.5)
) +
  facet_grid(~variable, scales="free") +
  scale_x_discrete(breaks=c("C", "Tr_DomWin", "Tr_DomLose", "Tr_ForWin", "Tr_ForLose"), labels=c("Control", "T1: UK \n sector wins", "T2: UK \n sector loses", "T3: German \n sector wins", "T4: German \n sector loses")) +
  geom_hline(data=meandf, aes(yintercept=mean), linetype=c("dashed","dashed"), size=.5,col=c(cols[2],cols[3])) +
  theme_bw() +
  theme(text = element_text(size=11),
        axis.text.x = element_text(size=9),
        legend.position="none"
  )
dev.off()

# (3) HEARD ABOUT PARIS AGREEMENT

f$paris <- ifelse(f$heardparisaccord<3,1,ifelse(f$heardparisaccord==3,2,NA))

# Restructure data for better plotting options
df.plot <- f[is.na(f$paris)==FALSE,]
df.plot$paris <- factor(df.plot$paris, levels=c(1,2), labels=c("Heard of Paris", "Not heard of Paris"))


# Beliefs for: Heard about Paris
group_by(df.plot[df.plot$paris=="Heard of Paris",], treat_vign) %>%
  summarise(
    count = n(),
    mean = mean(vign_beliefUKwins, na.rm = TRUE),
    sd = sd(vign_beliefUKwins, na.rm = TRUE)
  )

# Beliefs for: Has not heard about Paris
group_by(df.plot[df.plot$paris=="Not heard of Paris",], treat_vign) %>%
  summarise(
    count = n(),
    mean = mean(vign_beliefUKwins, na.rm = TRUE),
    sd = sd(vign_beliefUKwins, na.rm = TRUE)
  )

# Compare respondents across groups
# Multiple comparison p-value correction
pvallist <- NA
for (i in 1:5) {
  pvallist[i] <- t.test(vign_beliefUKwins ~ paris, data=df.plot[as.numeric(df.plot$treat_vign)==i ,])$p.value
}
p.adjust(pvallist,method="BH")



# ================================================================
# Appendix A7.8: Education
# ================================================================


f$edu <- ifelse(f$education<4,1,ifelse((f$education>=4 & f$education!=7),2,NA))
subdf <- f[is.na(f$edu)==FALSE,]
subdf$edu <- factor(subdf$edu, levels=c(1,2), labels=c("No university education", "University education"))


# Beliefs for no universiy education
h1 <- group_by(subdf[subdf$edu=="No university education",], treat_vign) %>%
  summarise(
    count = n(),
    mean = mean(vign_beliefUKwins, na.rm = TRUE),
    sd = sd(vign_beliefUKwins, na.rm = TRUE)
  )

# ANOVA: No university education
res.aov <- aov(vign_beliefUKwins ~ treat_vign, data = subdf[subdf$edu=="No university education",])
summary(res.aov)
TukeyHSD(res.aov)

# Beliefs for University education
h2 <- group_by(subdf[subdf$edu=="University education",], treat_vign) %>%
  summarise(
    count = n(),
    mean = mean(vign_beliefUKwins, na.rm = TRUE),
    sd = sd(vign_beliefUKwins, na.rm = TRUE)
  )

# ANOVA: University education
res.aov <- aov(vign_beliefUKwins ~ treat_vign, data = subdf[subdf$edu=="University education",])
summary(res.aov)
TukeyHSD(res.aov)

# Compare subgroups
# Multiple comparison p-value correction
pvallist <- NA
for (i in 1:5) {
  pvallist[i] <- t.test(vign_beliefUKwins ~ edu, data=subdf[as.numeric(subdf$treat_vign)==i,])$p.value
}
p.adjust(pvallist,method="BH")


# ================================================================
# Figure A15: No university education vs university education
# ================================================================


# Restructure data for better plotting options
subdf.plot <- cbind(subdf[,c("treat_vign","edu","vign_beliefUKwins")])
colnames(subdf.plot)[2:3] <- c("variable", "outcome")


pdf(file="./belief_UK_win_edu.pdf", width=8, height=5)

cols <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")

# Construct data set for means
meandf <- data.frame(variable=c(levels(subdf.plot$variable)), mean=c(h1$mean[1],h2$mean[1]))

ggerrorplot(subdf.plot, x="treat_vign", y="outcome", color="variable",
            desc_stat="mean_ci", 
            error.plot="errorbar",
            add="mean",
            legend="none",
            ylab="DV (ordered)", xlab = "Experimental conditions",
            palette=c(cols[2],cols[3]),
            ylim=c(1.5,3.5)
) +
  facet_grid(~variable, scales="free") +
  scale_x_discrete(breaks=c("C", "Tr_DomWin", "Tr_DomLose", "Tr_ForWin", "Tr_ForLose"), labels=c("Control", "T1: UK \n sector wins", "T2: UK \n sector loses", "T3: German \n sector wins", "T4: German \n sector loses")) +
  geom_hline(data=meandf, aes(yintercept=mean), linetype=c("dashed","dashed"), size=.5,col=c(cols[2],cols[3])) +
  theme_bw() +
  theme(text = element_text(size=11),
        axis.text.x = element_text(size=9),
        legend.position="none"
  )
dev.off()


# ================================================================
# Appendix A7.9: Country performance on climate action
# ================================================================

# (1) UK PERFORMANCE 

f$perf <- ifelse(f$UK_ccperformance>2,1,ifelse(f$UK_ccperformance<3,2,NA))

# Restructure data for better plotting options
subdf <- f[is.na(f$perf)==FALSE,]
subdf$perf <- factor(subdf$perf, levels=c(1,2), labels=c("UK performs well", "UK performs poorly"))


# Beliefs for those thinking of UK as good performer
h1 <- group_by(subdf[subdf$perf=="UK performs well",], treat_vign) %>%
  summarise(
    count = n(),
    mean = mean(vign_beliefUKwins, na.rm = TRUE),
    sd = sd(vign_beliefUKwins, na.rm = TRUE)
  )

# Beliefs for those thinking of UK as bad performer
h2 <- group_by(subdf[subdf$perf=="UK performs poorly",], treat_vign) %>%
  summarise(
    count = n(),
    mean = mean(vign_beliefUKwins, na.rm = TRUE),
    sd = sd(vign_beliefUKwins, na.rm = TRUE)
  )


# Compare subgroups
# Multiple comparison p-value correction
pvallist <- NA
for (i in 1:5) {
  pvallist[i] <- t.test(vign_beliefUKwins ~ perf, data=subdf[as.numeric(subdf$treat_vign)==i,])$p.value
}
p.adjust(pvallist,method="BH")


# ================================================================
# Figure A16: UK performs well vs UK performs poorly
# ================================================================


# Restructure data for better plotting options
subdf.plot <- cbind(subdf[,c("treat_vign","perf","vign_beliefUKwins")])
colnames(subdf.plot)[2:3] <- c("variable", "outcome")


pdf(file="./belief_UK_win_perf.pdf", width=8, height=5)

cols <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")

# Construct data set for means
meandf <- data.frame(variable=c(levels(subdf.plot$variable)), mean=c(h1$mean[1],h2$mean[1]))

ggerrorplot(subdf.plot, x="treat_vign", y="outcome", color="variable",
            desc_stat="mean_ci", 
            error.plot="errorbar",
            add="mean",
            legend="none",
            ylab="DV (ordered)", xlab = "Experimental conditions",
            palette=c(cols[2],cols[3]),
            ylim=c(1.5,3.5)
) +
  facet_grid(~variable, scales="free") +
  scale_x_discrete(breaks=c("C", "Tr_DomWin", "Tr_DomLose", "Tr_ForWin", "Tr_ForLose"), labels=c("Control", "T1: UK \n sector wins", "T2: UK \n sector loses", "T3: German \n sector wins", "T4: German \n sector loses")) +
  geom_hline(data=meandf, aes(yintercept=mean), linetype=c("dashed","dashed"), size=.5,col=c(cols[2],cols[3])) +
  theme_bw() +
  theme(text = element_text(size=11),
        axis.text.x = element_text(size=9),
        legend.position="none"
  )
dev.off()


# ================================================================
# Figure A17: Developed countries perform well vs developed countries perform poorly
# ================================================================

f$dperf <- ifelse(f$devep_ccperformance<5 & f$devep_ccperformance>2,1,ifelse(f$devep_ccperformance<3,2,NA))

# Restructure data for better plotting options
subdf <- f[is.na(f$dperf)==FALSE,]
subdf$dperf <- factor(subdf$dperf, levels=c(1,2), labels=c("Developed countries perform well", "Developed countries perform poorly"))



# Beliefs for those thinking of develeled countries as good performers
h1 <- group_by(subdf[subdf$dperf=="Developed countries perform well",], treat_vign) %>%
  summarise(
    count = n(),
    mean = mean(vign_beliefUKwins, na.rm = TRUE),
    sd = sd(vign_beliefUKwins, na.rm = TRUE)
  )

# Beliefs for those thinking of developed countries as bad performers
h2 <- group_by(subdf[subdf$dperf=="Developed countries perform poorly",], treat_vign) %>%
  summarise(
    count = n(),
    mean = mean(vign_beliefUKwins, na.rm = TRUE),
    sd = sd(vign_beliefUKwins, na.rm = TRUE)
  )


# Compare subgroups
# Multiple comparison p-value correction
pvallist <- NA
for (i in 1:5) {
  pvallist[i] <- t.test(vign_beliefUKwins ~ dperf, data=subdf[as.numeric(subdf$treat_vign)==i,])$p.value
}
p.adjust(pvallist,method="BH")


# ================================================================
# Figure A17: Developed countries perform well vs poorly
# ================================================================


# Restructure data for better plotting options
subdf.plot <- cbind(subdf[,c("treat_vign","dperf","vign_beliefUKwins")])
colnames(subdf.plot)[2:3] <- c("variable", "outcome")


pdf(file="./belief_UK_win_dperf.pdf", width=8, height=5)

cols <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")

# Construct data set for means
meandf <- data.frame(variable=c(levels(subdf.plot$variable)), mean=c(h1$mean[1],h2$mean[1]))

ggerrorplot(subdf.plot, x="treat_vign", y="outcome", color="variable",
            desc_stat="mean_ci", 
            error.plot="errorbar",
            add="mean",
            legend="none",
            ylab="DV (ordered)", xlab = "Experimental conditions",
            palette=c(cols[2],cols[3]),
            ylim=c(1.5,3.5)
) +
  facet_grid(~variable, scales="free") +
  scale_x_discrete(breaks=c("C", "Tr_DomWin", "Tr_DomLose", "Tr_ForWin", "Tr_ForLose"), labels=c("Control", "T1: UK \n sector wins", "T2: UK \n sector loses", "T3: German \n sector wins", "T4: German \n sector loses")) +
  geom_hline(data=meandf, aes(yintercept=mean), linetype=c("dashed","dashed"), size=.5,col=c(cols[2],cols[3])) +
  theme_bw() +
  theme(text = element_text(size=11),
        axis.text.x = element_text(size=9),
        legend.position="none"
  )
dev.off()


# ================================================================
#                           END OF FILE
# ================================================================
